// Copyright (C) 2011  Andrew H. Chan
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

#include "buneman.h"
#include "util.h"
//#include "approx.h"
#include <glpk.h>
#include <float.h>

using namespace std;

std::size_t hash_value(boost::tuple<uint,uint,int> const& t);

Buneman::Buneman(Input_data const& in_data) : in_data(in_data), tot_new_num_sites(0) { }

int Buneman::execute() {
    preprocess_data(in_data);

    compute_buneman_property();
    create_segment_graphs();
    solve_segment_ILPs();

    compute_all_buneman_property();
    construct_level1_graph();

    construct_level2_graph();

    return 0;
}

int Buneman::preprocess_data(Input_data const& in_data) {
    cout << "---Pruning columns" << endl;

    row_data.clear();
    row_data.resize(NUM_SEGS);

    uint num_rows = in_data.processed_data.size();
    uint num_cols = in_data.processed_data[0].size();

    // remove duplicate rows for each segment
    // note that row_data is indexed by segment THEN haplotype
    // but seg_data is indexed by haplotype THEN segment
    for(uint i=0; i < NUM_SEGS; ++i) {
        boost::unordered_set<cvector> seg_set;
        for(uint j=0; j < num_rows; ++j) {
            if(seg_set.find(in_data.seg_data[j][i]) == seg_set.end()) {
                seg_set.insert(in_data.seg_data[j][i]);
                row_data[i].push_back(in_data.seg_data[j][i]);
            }
        }
    }

    ///////////////////////
    // prune columns
    //mark columns to be removed
    all_remove_col.clear();
    all_remove_col.resize(num_cols, false);
    int base_indx = 0;
    for(uint i=0; i < NUM_SEGS; ++i) {
        uint num_cols_removed = 0;
        uint num_sites = row_data[i][0].size();
        vector<bool> remove_col(num_sites, false);

        for(uint j=0; j< num_sites; ++j) {
            if(remove_col[j] == false) {
                bool contains_both = false;
                for(uint l=0; l < row_data[i].size(); ++l) {
                    if(row_data[i][l][j] != row_data[i][0][j]) {
                        contains_both = true;
                        break;
                    }
                }
                if(!contains_both && (num_cols_removed < num_sites-1)) {
                    // non-segregating columns were already removed in
                    //    processed_data, so no columns should be removed here
                    remove_col[j] = true;
                    all_remove_col[base_indx + j] = true;
                    ++num_cols_removed;
                }
            }
        }

        for(uint j=0; j < num_sites; ++j) {
            if(remove_col[j] == false) {
                bool contains_mult_seg = false;
                uint contains_0 = 0; 
                uint contains_1 = 0;
                for(uint k=0; k < row_data[i].size(); ++k) {
                    if(contains_0 <= 2 && row_data[i][k][j] == 0) {
                        ++contains_0;
                    }
                    if(contains_1 <= 2 && row_data[i][k][j] == 1) {
                        ++contains_1;
                    }
                    if(contains_1 >= 2 && contains_0 >= 2) {
                        contains_mult_seg = true;
                        break;
                    }
                }
                if(!contains_mult_seg && (num_cols_removed < num_sites-1)) {
                    remove_col[j] = true;
                    all_remove_col[base_indx + j] = true;
                    ++num_cols_removed;
                }
            }
        } // Done removing columns with a single bit flip 

        for(uint j=0; j < num_sites-1; ++j) {
            if(remove_col[j] == false) {
                for(uint k=j+1; k < num_sites; ++k) {
                    if(remove_col[k] == false) {
                        bool col_iden = true;
                        for(uint l=0; l < row_data[i].size(); ++l) {
                            if(row_data[i][l][j] != row_data[i][l][k]) {
                                col_iden = false;
                                break;
                            }
                        }
                        if(col_iden == true && (num_cols_removed < num_sites-1)) {
                            remove_col[k] = true; // mark column to be deleted
                            all_remove_col[base_indx + k] = true;
                            ++num_cols_removed;
                        }
                    }
                }
            }
        }

        for(uint j=0; j < num_sites-1; ++j) {
            if(!remove_col[j]) {
                for(uint k=j+1; k < num_sites; ++k) {
                    bool complement = true;
                    if(!remove_col[k]){
                        for(uint l=0; l < row_data[i].size(); ++l) {
                            if(row_data[i][l][j] == row_data[i][l][k]) {
                                complement = false; 
                                break;
                            }
                        }
                        if(complement && (num_cols_removed < num_sites-1)) {
                            remove_col[k] = true;
                            all_remove_col[base_indx + k] = true;
                            ++num_cols_removed;
                        }
                    }
                }
            }
        }
        assert(num_cols_removed < num_sites);

        // After marking columns to be removed 
        uint new_num_sites = 0;
        for(uint ii = 0; ii < remove_col.size(); ++ii) {
            if(!remove_col[ii]) {
                ++new_num_sites;
            }
        }

        assert(new_num_sites > 0);
        for(uint m=0; m < row_data[i].size(); ++m) {
            uint pp = 0;
            for(uint n=0; n < num_sites; ++n) {
                if(!remove_col[n]) {
                    row_data[i][m][pp] = row_data[i][m][n];
                    assert(pp <= n);
                    ++pp;
                }
            }
            row_data[i][m].resize(new_num_sites);
        }
        base_indx += num_sites;
    }
    // done processing columns

    // get rid of duplicate segment configurations again in case removing
    //    the columns above has made some rows identical.
    for(uint i=0; i < NUM_SEGS; ++i) {
        boost::unordered_set<cvector> seg_set;
        uint index_counter = 0;
        for(uint j=0; j < row_data[i].size(); ++j) {
            if(seg_set.find(row_data[i][j]) == seg_set.end()) {
                seg_set.insert(row_data[i][j]);
                row_data[i][index_counter] = row_data[i][j];
                ++index_counter;
            }
        }
        row_data[i].resize(index_counter);
        assert(row_data[i].size() > 0);
    }

    tot_new_num_sites = 0;
    for(uint ii=0; ii<NUM_SEGS; ++ii) {
        assert(row_data[ii][0].size() > 0);
        tot_new_num_sites += row_data[ii][0].size();
    }

    // construct map_to_data
    // maps sites in pruned data to sites in original data
    map_to_data.resize(tot_new_num_sites);
    uint pp = 0;
    for(uint ii=0; ii<num_cols; ++ii) {
        if(all_remove_col[ii] == false) {
            map_to_data[pp] = ii; 
            ++pp;
        }
    }
    assert(pp == tot_new_num_sites);

    // construct segment_start, segment_end and site_belong
    // Note: segment_start and segment_end contain indexes AFTER column pruning
    segment_start.resize(NUM_SEGS);
    segment_end.resize(NUM_SEGS);
    site_belong.resize(tot_new_num_sites);
    uint mm = 0;
    for(uint ii=0; ii<NUM_SEGS; ++ii) {
        // segment_start holds the index of the first site in the segment
        segment_start[ii] = mm;
        for(uint jj=0; jj<row_data[ii][0].size(); ++jj) {
            // site_belong maps a site to which segment it's in
            site_belong[mm] = ii; 
            ++mm;
        }
        // segment_end holds the index one beyond the last site in the segment
        segment_end[ii] = mm;
    }
    assert(mm == tot_new_num_sites);

    return 0;
}

int Buneman::compute_buneman_property() {
    bprop.clear();
    bprop.resize(NUM_SEGS);

    cout << "---Computing Buneman property for each segment" << endl;
    for(uint i=0; i < NUM_SEGS; ++i) {
        uint num_sites = row_data[i][0].size();
        assert(num_sites > 0);
        bprop[i].resize(num_sites * (num_sites-1) / 2, 0);
    }

    for(uint i=0; i < NUM_SEGS; ++i) {
        uint num_sites = row_data[i][0].size();
        for(uint j=0; j < num_sites - 1; ++j) {
            for(uint k=j+1; k < num_sites; ++k) {
                // number of haps in segment i
                uint row_data_size = row_data[i].size();
                // check if buneman property is satisfied for sites (j,k)
                // where j < k
                // 00
                bprop[i][j * (num_sites-2) - j*(j-1) / 2 + k - 1] = 0;
                for(uint l=0; l < row_data_size; ++l) {
                    if(row_data[i][l][j] == 0 && row_data[i][l][k] == 0) {
                        bprop[i][j * (num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x8;
                        break;
                    }
                }
                // 01
                for(uint l=0; l < row_data_size; ++l) {
                    if(row_data[i][l][j] == 0 && row_data[i][l][k] == 1) {
                        bprop[i][j * (num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x4;
                        break;
                    }
                }
                // 10
                for(uint l=0; l < row_data_size; ++l) {
                    if(row_data[i][l][j] == 1 && row_data[i][l][k] == 0) {
                        bprop[i][j * (num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x2;
                        break;
                    }
                }
                // 11
                for(uint l=0; l < row_data_size; ++l) {
                    if(row_data[i][l][j] == 1 && row_data[i][l][k] == 1) {
                        bprop[i][j * (num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x1;
                        break;
                    }
                }
            }
        } // end iterating over sites
    } // end iterating over segments

    cout << "Done computing Buneman property for each segment" << endl;
    return 0;
}

int Buneman::create_segment_graphs() {
    cout << "---Creating Buneman graph for each segment" << endl;
    // create buneman graphs
    buneman_graphs.resize(NUM_SEGS);
    for(uint i=0; i < NUM_SEGS; ++i) {
        num_edges = 0;
        buneman_graphs[i] = Segment_graph(); 
        Segment_graph &cur_graph = buneman_graphs[i];

        boost::unordered_set<cvector> v_set;
        // add one vertex to vertex set and explore from there
        v_set.insert(row_data[i][0]);
        cur_graph.add_vertex(row_data[i][0]);
        buneman_neighbor(i, v_set, row_data[i][0], cur_graph);
    }
    // add terminal vertices for each buneman graph
    for(uint i=0; i<NUM_SEGS; ++i) {
        for(uint j=0; j<row_data[i].size(); ++j) {
            buneman_graphs[i].graph_data.terminal_vertices.push_back(row_data[i][j]);
            // make sure every terminal vertex is unique
            assert(buneman_graphs[i].graph_data.terminal_set.find(row_data[i][j]) == buneman_graphs[i].graph_data.terminal_set.end());
            buneman_graphs[i].graph_data.terminal_set.insert(row_data[i][j]);
        }
        assert(buneman_graphs[i].num_vertices() >= buneman_graphs[i].graph_data.terminal_vertices.size());
    }
    return 0;
}

int Buneman::buneman_neighbor(uint seg_id, boost::unordered_set<cvector>& v_set, cvector const& vertex, Segment_graph& cur_graph) {
    uint num_sites = row_data[seg_id][0].size();
    for(uint i=0; i < num_sites; ++i) {
        cvector v = vertex;
        // flip the bit
        v[i] = v[i] ? 0 : 1;
        if(v_set.find(v) != v_set.end()) {
            if(cur_graph.vtoid[v] > cur_graph.vtoid[vertex]) {
                cur_graph.add_edge(vertex, v);
                cur_graph.add_edge(v, vertex);
            }
            continue;
        } else {
            // check if v is buneman
            bool v_is_buneman = true;
            for(uint j=0; j < num_sites - 1; ++j) {
                for(uint k=j+1; k < num_sites; ++k) {
                    char current_vertex = 0;
                    if(v[j] == 0 && v[k] == 0) {
                        current_vertex |= 0x8;
                    } else if(v[j] == 0 && v[k] == 1){
                        current_vertex |= 0x4;
                    } else if(v[j] == 1 && v[k] == 0){
                        current_vertex |= 0x2;
                    } else  {
                        assert(v[j] == 1 && v[k] == 1);
                        current_vertex |= 0x1;
                    }
                    if((current_vertex | bprop[seg_id][j * (num_sites-2) - j*(j-1) / 2 + k - 1]) != bprop[seg_id][j * (num_sites-2) - j*(j-1) / 2 + k - 1]) {
                        v_is_buneman = false;
                        break;
                    }
                } // end iterating over second site in pair
                if(!v_is_buneman) {
                    break;
                }
            } // end iterating over first site in pair
            if(num_sites == 1) {
                v_is_buneman = false;
                for(uint hapnum=0; hapnum < row_data[i].size(); ++hapnum) {
                    if(row_data[seg_id][hapnum][0] == v[0]) {
                        v_is_buneman = true;
                        break;
                    }
                }
            }
            if(v_is_buneman) {
                v_set.insert(v);
                pair<uint, bool> add_ret = cur_graph.add_vertex(v);
                assert(add_ret.second);
                cur_graph.add_edge(vertex, v);
                cur_graph.add_edge(v, vertex);
                buneman_neighbor(seg_id, v_set, v, cur_graph);
            }
        }
    } // end flipping each site
    return 0;
}

int Buneman::solve_segment_ILPs() {
    cout << "---Creating ILP for Buneman graph for each segment" << endl;
    steiner_vertices.clear();
    steiner_vertices.resize(NUM_SEGS);
    // for each buneman graph, solve the steiner tree problem
    // create MIP instance for each one
    for(uint i=0; i<NUM_SEGS; ++i) {
        Segment_graph &cur_graph = buneman_graphs[i];
        string lp_name = "buneman_" + boost::lexical_cast<string>(i);
        if(cur_graph.graph_data.terminal_vertices.size() > 1) {
            glp_prob *lp = glp_create_prob();
            glp_set_prob_name(lp, lp_name.c_str());

            /////////////////////////////////////////////////////
            // add variables of lp
            // add s variables for every edge
            boost::unordered_map<pair<uint,uint>, int> map_edge_svar;
            vector<boost::unordered_map<pair<uint,uint>, int> > map_edge_fvar;
            //map_edge_fvar.resize(cur_graph.graph_data.terminal_vertices.size());

            for(Segment_graph::adj_list::const_iterator viter = cur_graph.adj.begin(); viter != cur_graph.adj.end(); ++viter) {
                for(Segment_graph::edge_list::const_iterator eiter = viter->second.first.begin(); eiter != viter->second.first.end(); ++eiter) {
                    int s_id = glp_add_cols(lp, 1);
                    glp_set_col_kind(lp, s_id, GLP_BV); // binary
                    map_edge_svar[pair<uint,uint>(viter->first,eiter->first)] = s_id;
                }
            }
            uint root_id = cur_graph.vtoid.at(cur_graph.graph_data.terminal_vertices[0]);

            for(uint term_vert=1; term_vert < cur_graph.graph_data.terminal_vertices.size(); ++term_vert) {
                uint terminal_vertex_id = cur_graph.vtoid.at(cur_graph.graph_data.terminal_vertices.at(term_vert));
                boost::unordered_map<pair<uint, uint>, int> map_edge_flow;
                // for every terminal, create its flow variables 
                for(Segment_graph::adj_list::const_iterator viter=cur_graph.adj.begin(); viter!=cur_graph.adj.end(); ++viter) {
                    for(Segment_graph::edge_list::const_iterator eiter = viter->second.first.begin(); eiter != viter->second.first.end(); ++eiter) {
                        int s_id = map_edge_svar.at(pair<uint,uint>(viter->first, eiter->first));
                        int flow_id = glp_add_cols(lp, 1);
                        //map_edge_fvar[term_vert][pair<uint,uint>(v_id, *eiter)] = flow_id;
                        map_edge_flow[pair<uint, uint>(viter->first, eiter->first)] = flow_id;
                        // add column constraint: 0 <= f_uv
                        //glp_set_col_bnds(lp, flow_id, GLP_LO, 0.0, 0.0);

                        // make flow binary
                        glp_set_col_kind(lp, flow_id, GLP_BV);

                        // add row constraint: f_uv <= s_uv
                        //                     0 <= s_uv - f_uv
                        int row_id = glp_add_rows(lp, 1);
                        glp_set_row_bnds(lp, row_id, GLP_LO, 0.0, 0.0);
                        int ind[3] = {0, s_id, flow_id};
                        double val[3] = {0.0, 1.0, -1.0};
                        glp_set_mat_row(lp, row_id, 2, ind, val);
                    }
                } // done adding s variables and flow variables for terminal t

                // add flow conservation constraints for flow t 
                for(Segment_graph::adj_list::const_iterator viter=cur_graph.adj.begin(); viter != cur_graph.adj.end(); ++viter) {
                    // don't need flow conservation for root
                    if(viter->first == root_id) continue;
                    // don't need flow conservation for terminal t 
                    //if(cur_graph.terminal_set.find(cur_graph.idtov[v_id]) != cur_graph.terminal_set.end()) continue;
                    if(viter->first == terminal_vertex_id) continue;

                    // for v_id, create flow constraints for vertex
                    int conserve_id = glp_add_rows(lp, 1);
                    glp_set_row_bnds(lp, conserve_id, GLP_FX, 0.0, 0.0);

                    int len = viter->second.first.size() + viter->second.second.size();
                    int ind[len+1]; 
                    double val[len+1];

                    int index_counter = 1;
                    // iterate over outgoing edges
                    for(Segment_graph::edge_list::const_iterator eiter = viter->second.first.begin(); eiter != viter->second.first.end(); ++eiter) {
                        uint fromedge = viter->first;
                        uint toedge = eiter->first; 
                        // outgoing edge
                        assert(map_edge_flow.find(pair<uint,uint>(fromedge,toedge)) != map_edge_flow.end());
                        ind[index_counter] = map_edge_flow[pair<uint,uint>(fromedge, toedge)];
                        val[index_counter] = 1.0;
                        ++index_counter;
                    }
                    // iterate over incoming edges
                    for(Segment_graph::edge_list::const_iterator eiter = viter->second.second.begin(); eiter != viter->second.second.end(); ++eiter) {
                        uint fromedge = eiter->first; 
                        uint toedge = viter->first;
                        // incoming edge
                        assert(map_edge_flow.find(pair<uint,uint>(fromedge,toedge)) != map_edge_flow.end());
                        ind[index_counter] = map_edge_flow[pair<uint,uint>(fromedge, toedge)];
                        val[index_counter] = -1.0;
                        ++index_counter;
                    }
                    assert(index_counter == len+1);

                    glp_set_mat_row(lp, conserve_id, len, ind, val);

                } // done adding flow constraints for flow for term t

                // add constraint to prevent flow from leaving terminal t
                int term_id = glp_add_rows(lp, 1);
                glp_set_row_bnds(lp, term_id, GLP_FX, 0.0, 0.0);

                int len = cur_graph.adj.at(terminal_vertex_id).first.size();

                int ind[len+1];
                double val[len+1];

                int index_counter = 1;
                // iterate over outgoing edges
                for(Segment_graph::edge_list::const_iterator eiter = cur_graph.adj.at(terminal_vertex_id).first.begin(); eiter != cur_graph.adj.at(terminal_vertex_id).first.end(); ++eiter) {
                    uint fromedge = terminal_vertex_id; 
                    uint toedge = eiter->first; 
                    // outgoing edge
                    assert(map_edge_flow.find(pair<uint,uint>(fromedge,toedge)) != map_edge_flow.end());
                    ind[index_counter] = map_edge_flow[pair<uint,uint>(fromedge, toedge)];
                    val[index_counter] = 1.0;
                    ++index_counter;
                }
                assert(index_counter == len+1);

                glp_set_mat_row(lp, term_id, len, ind, val);
                // done adding term outflow constraints for flow t 

                // add constraint so that terminal receives 1 unit of flow 
                int term_in_id = glp_add_rows(lp, 1);
                glp_set_row_bnds(lp, term_in_id, GLP_FX, 1.0, 1.0);

                int len_in = cur_graph.adj[terminal_vertex_id].second.size();

                int ind_in[len_in+1]; 
                double val_in[len_in+1];

                index_counter = 1;
                // iterate over incoming edges
                for(Segment_graph::edge_list::const_iterator eiter = cur_graph.adj.at(terminal_vertex_id).second.begin(); eiter != cur_graph.adj.at(terminal_vertex_id).second.end(); ++eiter) {
                    uint fromedge = eiter->first; 
                    uint toedge = terminal_vertex_id; 
                    // incoming edge
                    assert(map_edge_flow.find(pair<uint,uint>(fromedge,toedge)) != map_edge_flow.end());
                    ind_in[index_counter] = map_edge_flow[pair<uint,uint>(fromedge, toedge)];
                    val_in[index_counter] = 1.0;
                    ++index_counter;
                }
                assert(index_counter == len_in+1);

                glp_set_mat_row(lp, term_in_id, len_in, ind_in, val_in);
                // done adding term inflow constraints 

            } // done processing the flow for terminal term_vert

            // set objective function
            for(boost::unordered_map<pair<uint,uint>, int>::const_iterator siter = map_edge_svar.begin(); siter != map_edge_svar.end(); ++siter) {
                glp_set_obj_coef(lp, siter->second, 1.0);
            }

            // solve mip
            glp_iocp parm;
            glp_init_iocp(&parm); // init parameters struct 
            parm.presolve = GLP_ON; // turn on presolver
            parm.msg_lev = GLP_MSG_ERR; // only output error messages 

            int ret = glp_intopt(lp, &parm);
            assert(ret == 0);

            // display solution, s vars
            //cout << "svars" << endl;
            //for(boost::unordered_map<pair<uint,uint>, int>::const_iterator siter = map_edge_svar.begin(); siter != map_edge_svar.end(); ++siter) {
            //    double sol_value = glp_mip_col_val(lp, siter->second);
            //    if(sol_value != 0) {
            //        cout << "edge: ";
            //        cout << cur_graph.idtov[siter->first.first] << " " << cur_graph.idtov[siter->first.second] << endl;
            //        cout << "val:  " << sol_value << endl;
            //    }
            //}

            // display solution, f vars
            //cout << "flow vars" << endl;
            //cout << "root: " << cur_graph.terminal_vertices[0] << endl;
            //for(uint tt=1; tt<cur_graph.terminal_vertices.size(); ++tt) {
            //    cout << "flow for terminal " << tt << endl;
            //    cout << "terminal: " << cur_graph.terminal_vertices[tt] << endl;
            //    for(boost::unordered_map<pair<uint,uint>, int>::const_iterator fiter = map_edge_fvar[tt].begin(); fiter != map_edge_fvar[tt].end(); ++fiter) {
            //        double sol_value = glp_mip_col_val(lp, fiter->second);
            //        if(sol_value != 0) {
            //            cout << "edge: ";
            //            cout << cur_graph.idtov[fiter->first.first] << " " << cur_graph.idtov[fiter->first.second] << endl;
            //            cout << "val:  " << sol_value << endl;
            //        }
            //    }
            //}

            //int obj_value = glp_mip_obj_val(lp);
            //cout << "obj_value: " << obj_value << endl;
            // this is a bad assertion since steiner tree might use
            // intermediate vertices
            //assert(obj_value == static_cast<int>(cur_graph.terminal_vertices.size()) - 1);

            // record steiner vertices
            steiner_vertices[i].push_back(cur_graph.idtov.at(root_id));
            for(boost::unordered_map<pair<uint,uint>, int>::const_iterator siter = map_edge_svar.begin(); siter != map_edge_svar.end(); ++siter) {
                double sol_value = glp_mip_col_val(lp, siter->second);
                if(sol_value != 0) {
                    steiner_vertices[i].push_back(cur_graph.idtov.at(siter->first.second));
                }
            }
            assert(steiner_vertices[i].size() >= cur_graph.graph_data.terminal_vertices.size());

            // delete lp for segment
            glp_delete_prob(lp);
        } else {
            // this is to handle the case where there is
            //    only one terminal vertex
            steiner_vertices[i].push_back(cur_graph.graph_data.terminal_vertices[0]);
        }
        cout << "Found Steiner tree for segment " << i << endl;
    }
    // done with ILP for steiner trees for individual buneman graphs

    for(uint ii=0; ii < NUM_SEGS; ++ii) {
        cout << "size of steiner tree " << ii << ": " << steiner_vertices[ii].size() << ", num terminal vertices: " <<  buneman_graphs[ii].graph_data.terminal_vertices.size() << endl;
    }

    return 0;
}

int Buneman::compute_all_buneman_property() {
    ///////////////////////////////////////////////////////////////////////
    // preprocess processed_data in preparation of creating the main graph

    uint num_cols = in_data.processed_data[0].size();
    assert(num_cols > 0);
    uint num_haps = in_data.processed_data.size();

    allbprop.resize(tot_new_num_sites * (tot_new_num_sites-1) / 2, 0);
    allbprop_count.resize(tot_new_num_sites * (tot_new_num_sites-1) / 2, 0);

    for(uint j=0; j < tot_new_num_sites-1; ++j) {
        int index_j = map_to_data[j]; 
        for(uint k=j+1; k < tot_new_num_sites; ++k) {
            int index_k = map_to_data[k]; 
            int bprop_count = 0;
            // check if buneman property is satisfied for sites (j,k)
            // where j < k
            // 00
            for(uint l=0; l < num_haps; ++l) {
                if(in_data.processed_data[l][index_j] == 0 && in_data.processed_data[l][index_k] == 0) {
                    allbprop[j * (tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x8;
                    ++bprop_count;
                    break;
                }
            }
            // 01
            for(uint l=0; l < num_haps; ++l) {
                if(in_data.processed_data[l][index_j] == 0 && in_data.processed_data[l][index_k] == 1) {
                    allbprop[j * (tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x4;
                    ++bprop_count;
                    break;
                }
            }
            // 10
            for(uint l=0; l < num_haps; ++l) {
                if(in_data.processed_data[l][index_j] == 1 && in_data.processed_data[l][index_k] == 0) {
                    allbprop[j * (tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x2;
                    ++bprop_count;
                    break;
                }
            }
            // 11
            for(uint l=0; l < num_haps; ++l) {
                if(in_data.processed_data[l][index_j] == 1 && in_data.processed_data[l][index_k] == 1) {
                    allbprop[j * (tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] |= 0x1;
                    ++bprop_count;
                    break;
                }
            }
            // update count
            allbprop_count[j * (tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] = bprop_count;
            assert(bprop_count >= 0);
            assert(bprop_count <= 4);
        }
    } // end iterating over sites
    cout << "allbprop_count: " << allbprop_count << endl;
    return 0;
}

int Buneman::construct_level1_graph() {
    // Construct level 1 Buneman graph
    cout << "---Constructing level 1 graph" << endl;

    uint tot_new_num_sites1 = 0;
    for(uint i=0; i<NUM_SEGS; ++i) {
        tot_new_num_sites1 += steiner_vertices[i][0].size();
    }

    uint tot_new_num_sites2 = 0;
    for(uint i=0; i<NUM_SEGS; ++i) {
        tot_new_num_sites2 += row_data[i][0].size();
    }
    assert(tot_new_num_sites == tot_new_num_sites1);
    assert(tot_new_num_sites1 == tot_new_num_sites2);

    cout << "tot num sites:  " << tot_new_num_sites << endl;
    cout << "orig num sites: " << in_data.processed_data[0].size() << endl;

    vertices_not_kept = 0;
    num_vertices_added = 0;
    cvector tmp_part_vertex;

    ////////////////////////////////
    // Add vertices to main graph
    // iterate over possible first segments
    // add_level1_vertices recursively iterates over other segments
    add_level1_vertices(level1_graph, 0, tmp_part_vertex);
    cout << "num vertices added: " << num_vertices_added << endl;
    cout << "vertices not kept: " << vertices_not_kept << endl;

    // create list of terminal vertices 
    assert(in_data.processed_data[0].size() == all_remove_col.size());
    //boost::unordered_set<cvector> terminal_set;
    //  generate list of corresponding terminal vertices
    for(uint ii=0; ii < in_data.processed_data.size(); ++ii) {
        cvector tmp_hap(tot_new_num_sites, 0);
        uint mm = 0;
        for(uint jj=0; jj < in_data.processed_data[ii].size(); ++jj ) {
            if(all_remove_col[jj]==false) {
                tmp_hap[mm] = in_data.processed_data[ii][jj];
                ++mm;
            }
        }
        assert(mm == tot_new_num_sites);

        // after removing columns, some rows may become identical
        // make sure that every terminal vertex is unique 
        if(level1_graph.graph_data.terminal_set.find(tmp_hap) == level1_graph.graph_data.terminal_set.end()) {
            level1_graph.graph_data.terminal_set.insert(tmp_hap);
            level1_graph.graph_data.terminal_vertices.push_back(tmp_hap);
        }
    }
    //assert(main_graph.terminal_vertices.size() == in_data.processed_data.size());
    cout << "number of terminal vertices: " << level1_graph.graph_data.terminal_vertices.size() << endl;

    // check that terminal vertices are in the main graph
    for(uint ii=0; ii<level1_graph.graph_data.terminal_vertices.size(); ++ii) {
        assert(level1_graph.vtoid.find(level1_graph.graph_data.terminal_vertices[ii]) != level1_graph.vtoid.end());
    }

    cout << "root vertex: " << level1_graph.graph_data.terminal_vertices[0] << endl;
    cout << "---" << endl;
    cout << "number of terminal vertices: " << level1_graph.graph_data.terminal_vertices.size() << endl;
    cout << "terminal vertices of main graph" << endl;
    for(uint ii=0; ii<level1_graph.graph_data.terminal_vertices.size(); ++ii) {
        cout << level1_graph.graph_data.terminal_vertices[ii] << endl;
    }

    cout << "num vertices in level 1 graph: " << level1_graph.num_vertices() << endl;

    ////////////////////////////////////////
    // add edges to graph
    add_level1_edges(level1_graph);

    cout << "num edges in level 1 graph: " << level1_graph.num_edges() << endl;
    assert(level1_graph.check_multiedge() == 0);

    // check for connectedness
    boost::unordered_set<cvector> v_set;
    int num_vert_reached = 0;
    dfs_search(level1_graph, level1_graph.graph_data.terminal_vertices[0], v_set, num_vert_reached);
    cout << "num vert reached by dfs: " << num_vert_reached << endl;
    assert(num_vert_reached == static_cast<int>(level1_graph.num_vertices()));

    return 0;
}

int Buneman::add_level1_vertices(Level1_graph& graph, int segment, cvector const& part_vertex) {
    for(uint vert_counter=0; vert_counter < row_data[segment].size(); ++vert_counter) {
        cvector tmp_part_vertex = append_vectors(part_vertex, row_data[segment][vert_counter]);
        bool keep_going = true;
        // check for pruning 4th binary pair
        for(int j=0; j < segment_start[segment]; ++j) {
            for(int k=segment_start[segment]; k < segment_end[segment]; ++k) {
                assert(site_belong[j] != site_belong[k]);
                if(allbprop_count[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] == 3) {
                    char current_pair = tmp_part_vertex[j]?(tmp_part_vertex[k]?0x1:0x2):(tmp_part_vertex[k]?0x4:0x8);
                    if((current_pair | allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) != allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) {
                        keep_going = false;
                        break;
                    }
                }
            } // end iterating over k
            if(!keep_going) {
                break;
            }
        } // end iterating over j

        if(keep_going) {
            if(segment == static_cast<int>(NUM_SEGS)-1) { // base case
                pair<uint, bool> add_ret = graph.add_vertex(tmp_part_vertex);
                assert(add_ret.second);
                ++num_vertices_added;
            } else {
                add_level1_vertices(graph, segment+1, tmp_part_vertex);
            }
        } else {
            // compute number of vertices saved
            int num_vertices_saved = 1;
            for(uint qq=segment+1; qq < NUM_SEGS; ++qq) {
                num_vertices_saved *= row_data[qq].size();
            }
            vertices_not_kept += num_vertices_saved;
        }
    }
    return 0;
}

int Buneman::add_level1_edges(Level1_graph& graph) {
    // go through each vertex in the graph, flip each bit, and see if
    //    the new vertex exists
    // if it does exist, then create an edge
    for(Level1_graph::v_to_id::const_iterator viter = graph.vtoid.begin(); viter != graph.vtoid.end(); ++viter) {
        // flip bits one by one
        for(uint jj=0; jj<tot_new_num_sites; ++jj) {
            cvector tmp_v = viter->first; 
            tmp_v[jj] = viter->first[jj]?0:1;
            if(graph.vtoid.find(tmp_v) != graph.vtoid.end()) {
                graph.add_edge(viter->first, tmp_v);
            }
        }
    }
    return 0;
}

int Buneman::construct_level2_graph() {
    cout << "---Constructing level 2 graph" << endl;

    // Set costs for mutation and reassortment
    level2_graph.graph_data.edge_costs[0] = 1.0; // mutation
    level2_graph.graph_data.edge_costs[1] = 1.0; // reassortment

    // Construct graph
    vertices_not_kept2 = 0;
    dcvector tmp_part_vertex;
    // 4th parameter is heuristics
    add_level2_vertices(level2_graph, 0, tmp_part_vertex, 0);
    cout << "num vertices in level 2 graph: " << level2_graph.num_vertices() << endl;
    add_level2_edges(level2_graph);

    // add terminal vertices
    for(vector<Level1_graph::vertex_type>::const_iterator titer = level1_graph.graph_data.terminal_vertices.begin(); titer != level1_graph.graph_data.terminal_vertices.end(); ++titer) {
        level2_graph.graph_data.terminal_vertices.push_back(dcvector(*titer,*titer));
        level2_graph.graph_data.terminal_set.insert(dcvector(*titer,*titer));
        assert(level2_graph.vtoid.find(dcvector(*titer,*titer)) != level2_graph.vtoid.end());
    }
    // check for connectedness
    boost::unordered_set<dcvector> v_set;
    int num_vert_reached = 0;
    dfs_search(level2_graph, level2_graph.graph_data.terminal_vertices.at(0), v_set, num_vert_reached);
    cout << "num vert reached by dfs: " << num_vert_reached << endl;
    return 0;
}

    //create_big_graph();
    //cout << "num vertices in big graph: " << big_graph.num_vertices() << endl;
    //cout << "num edges in big graph: " << big_graph.num_edges() << endl;

    //big_graph.edge_costs[0] = 1.0; // mutation cost on left hap
    //big_graph.edge_costs[1] = 1.0; // mutation cost on right hap
    //big_graph.edge_costs[2] = 0.1; // reassortment cost on left hap
    //big_graph.edge_costs[3] = 0.1; // reassortment cost on right hap
    //big_graph.edge_costs[4] = 0.0; // special reassortment
    //big_graph.edge_costs[5] = 0.0; // special reassortment 

    ////create_reassort_ILP(big_graph);

    //cout << "approx steiner" << endl;
    //SteinerApprox steineralg;

    //Graph<pair<cvector, cvector> > modified_graph = big_graph;
    //modified_graph.edge_costs[6] = 0.0;
    //modified_graph.terminal_vertices.clear();
    //modified_graph.terminal_set.clear();
    //for(uint ii=0; ii < big_graph.terminal_vertices.size(); ++ii) {
    //    // create dummy terminals
    //    cvector dummy_hap;
    //    for(uint jj=0; jj < big_graph.terminal_vertices[ii].first.size(); ++jj) {
    //        if(big_graph.terminal_vertices[ii].first[jj] == 0){
    //            dummy_hap.push_back(2);
    //        } else {
    //            dummy_hap.push_back(3);
    //        }
    //    }
    //    pair<cvector,cvector> dummy_terminal;
    //    dummy_terminal = pair<cvector,cvector>(dummy_hap, dummy_hap);
    //    modified_graph.terminal_vertices.push_back(dummy_terminal);
    //    modified_graph.terminal_set.insert(dummy_terminal);
    //    pair<uint, bool> add_ret = modified_graph.add_unique_vertex(dummy_terminal);
    //    assert(add_ret.second);
    //    modified_graph.add_edge(big_graph.terminal_vertices[ii], dummy_terminal, 6, 0,0);
    //}
    //cout << "terminals in modified graph" << endl;
    //for(boost::unordered_set<pair<cvector, cvector> >::const_iterator viter = modified_graph.terminal_set.begin(); viter != modified_graph.terminal_set.end(); ++viter) {
    //    cout << viter->first << endl;
    //}

    //cout << "running steiner alg" << endl;
    //cout << "num vertices of modified graph: " << modified_graph.num_vertices() << endl;
    //cout << "num edges of modified graph: " << modified_graph.num_edges() <<  endl;
    //Graph<pair<cvector, cvector> > steinertree = steineralg.steiner_approx(modified_graph, big_graph.terminal_vertices[0], modified_graph.terminal_set, 2, modified_graph.terminal_set.size());
    //steineralg.shortest_path<pair<cvector,cvector> >(big_graph, big_graph.terminal_vertices[0], big_graph.terminal_vertices[1]);
    //cout << "num vertices in steiner tree: " << steinertree.num_vertices() << endl;
    //cout << "num edges in steiner tree: " << steinertree.num_edges() << endl;

    //for(Graph<pair<cvector,cvector> >::v_to_id::const_iterator viter = steinertree.vtoid.begin(); viter != steinertree.vtoid.end(); ++viter) {
    //    cout << viter->first << endl;
    //}
    //for(uint ii=0; ii < big_graph.terminal_vertices.size(); ++ii) {
    //    assert(steinertree.vtoid.find(big_graph.terminal_vertices[ii]) != steinertree.vtoid.end());
    //}
    //cout << "edges of steiner tree" << endl;
    //for(adj_list::const_iterator viter = steinertree.adj.begin(); viter != steinertree.adj.end(); ++viter) {
    //    for(edge_list::const_iterator eiter=viter->second.first.begin(); eiter!=viter->second.first.end(); ++eiter) {
    //        cout << steinertree.idtov[viter->first] << steinertree.idtov[eiter->get<0>()] << endl;
    //}

int Buneman::add_level2_vertices(Level2_graph& graph, int segment, dcvector const& part_vertex, int heuristics) {
    for(uint vert_counter1 = 0; vert_counter1 < row_data[segment].size(); ++vert_counter1) {
        dcvector tmp_vertex_first = dcvector(append_vectors(part_vertex.first, row_data[segment][vert_counter1]), part_vertex.second);
        for(uint vert_counter2 = 0; vert_counter2 < row_data[segment].size(); ++vert_counter2) {
            dcvector v = dcvector(tmp_vertex_first.first, append_vectors(part_vertex.second, row_data[segment][vert_counter2]));

            bool keep_going = true;

            ///////////////////////////////////////
            // check for pruning 4th binary pair

            // first coordinate
            for(int j=0; j < segment_start[segment]; ++j) {
                for(int k=segment_start[segment]; k < segment_end[segment]; ++k) {
                    assert(site_belong[j] != site_belong[k]);
                    if(allbprop_count[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] == 3) {
                        char current_pair = v.first[j]?(v.first[k]?0x1:0x2):(v.first[k]?0x4:0x8);
                        if((current_pair | allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) != allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) {
                            keep_going = false;
                            break;
                        }
                    }
                } // end iterating over k
                if(!keep_going) {
                    break;
                }
            } // end iterating over j
            // end checking for 4th binary pair in first coordinate

            // second coordinate
            if(keep_going) {
                for(int j=0; j < segment_start[segment]; ++j) {
                    for(int k=segment_start[segment]; k < segment_end[segment]; ++k) {
                        assert(site_belong[j] != site_belong[k]);
                        if(allbprop_count[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1] == 3) {
                            char current_pair = v.second[j]?(v.second[k]?0x1:0x2):(v.second[k]?0x4:0x8);
                            if((current_pair | allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) != allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) {
                                keep_going = false;
                                break;
                            }
                        }
                    } // end iterating over k
                    if(!keep_going) {
                        break;
                    }
                } // end iterating over j
            } // end checking for 4th binary pair in second coordinate

            ///////////////////////////////////////////////
            // The two filters below are heuristics
            // (they might prune away the optimal solutions)
            // check left (first) against right (second)
            if(heuristics == 1) {
                if(keep_going) {
                    for(int j=0; j < segment_start[segment]; ++j) {
                        for(int k=segment_start[segment]; k < segment_end[segment]; ++k) {
                            assert(site_belong[j] != site_belong[k]);
                            char current_pair = v.first[j]?(v.second[k]?0x1:0x2):(v.second[k]?0x4:0x8);
                            if((current_pair | allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) != allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) {
                                keep_going = false;
                                break;
                            }
                        } // end iterating over k
                        if(!keep_going) {
                            break;
                        }
                    } // end iterating over j
                }

                if(keep_going) {
                    // check right against left
                    for(int j=0; j < segment_start[segment]; ++j) {
                        for(int k=segment_start[segment]; k < segment_end[segment]; ++k) {
                            assert(site_belong[j] != site_belong[k]);
                            char current_pair = v.second[j]?(v.first[k]?0x1:0x2):(v.first[k]?0x4:0x8);
                            if((current_pair | allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) != allbprop[j*(tot_new_num_sites-2) - j*(j-1) / 2 + k - 1]) {
                                keep_going = false;
                                break;
                            }
                        } // end iterating over k
                        if(!keep_going) {
                            break;
                        }
                    } // end iterating over j
                }
            }
            // end heuristics
            /////////////////////////////////////

            if(keep_going) {
                if(segment == static_cast<int>(NUM_SEGS)-1) { // base case
                    pair<uint, bool> add_ret = level2_graph.add_vertex(v);
                    assert(add_ret.second);
                } else { // inductive case
                    add_level2_vertices(graph, segment+1, v, heuristics);
                }
            } else {
                // compute number of vertices saved
                int num_vertices_saved = 1;
                for(uint qq=segment+1; qq < NUM_SEGS; ++qq) {
                    num_vertices_saved *= row_data[qq].size();
                }
                // this could be an overestimate
                num_vertices_saved *= num_vertices_saved;
                vertices_not_kept2 += num_vertices_saved;
            }
        } // end vert_counter2
    } // end vert_counter1
    return 0;
}

int Buneman::add_level2_edges(Level2_graph& graph) {
    cout << "Adding edges to level 2 graph" << endl;
    cout << "   Adding mutation edges" << endl;
    uint num_mut_edges = 0;

    uint num_sites = graph.vtoid.begin()->first.first.size();

    for(Level2_graph::v_to_id::const_iterator iter = graph.vtoid.begin(); iter != graph.vtoid.end(); ++iter) {
        // mutate first hap
        for(uint mm=0; mm < num_sites; ++mm) {
            dcvector proposed = iter->first;
            proposed.first[mm] = 1 - proposed.first[mm];

            if(graph.vtoid.find(proposed) != graph.vtoid.end()) {
                uint vid1 = graph.vtoid.at(iter->first);
                uint vid2 = graph.vtoid.at(proposed);
                int weight = 1;
                graph.add_edge(vid1, vid2, Level2_edge_data(0, weight, pair<char,char>(0, mm)));
                ++num_mut_edges;
            }
        }

        for(uint mm=0; mm < num_sites; ++mm) {
            dcvector proposed = iter->first;
            proposed.second[mm] = 1 - proposed.second[mm];

            if(graph.vtoid.find(proposed) != graph.vtoid.end()) {
                uint vid1 = graph.vtoid.at(iter->first);
                uint vid2 = graph.vtoid.at(proposed);
                int weight = 1;
                graph.add_edge(vid1, vid2, Level2_edge_data(0, weight, pair<char,char>(1, mm)));
                ++num_mut_edges;
            }
        }
    }
    cout << "Number of mutation edges: " << num_mut_edges << endl;
    assert(graph.check_multiedge() == 0);

    // add reassortment edges
    cout << "   Adding reassortment edges" << endl;
    uint num_reassort_edges = 0;

    for(Level2_graph::v_to_id::const_iterator iter = graph.vtoid.begin(); iter != graph.vtoid.end(); ++iter) {
        vector<bool> mask(NUM_SEGS, false);
        while(true) {
            for(uint ii=0; ii<=mask.size(); ++ii) {
                if(mask[ii]) mask[ii] = false;
                else {
                    mask[ii] = true;
                    break;
                }
            }

            bool allones = true;
            for(uint ii=0; ii < mask.size(); ++ii) {
                if(mask[ii] == false) {
                    allones = false;
                    break;
                }
            }
            if(allones) break;

            // reassort first hap
            dcvector proposed1 = iter->first;
            for(uint kk=0; kk<NUM_SEGS; ++kk) {
                if(mask[kk]) {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj)
                        proposed1.first[jj] = iter->first.first[jj];
                } else {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj)
                        proposed1.first[jj] = iter->first.second[jj];
                }
            }


            char reassort_detail1 = 0;
            for(uint ll=0; ll < mask.size(); ++ll) {
                if(mask[ll]) {
                    reassort_detail1 |= (0x1 << ll);
                }
            }
            assert(reassort_detail1 != 0);

            // reassort second hap
            dcvector proposed2 = iter->first;
            for(uint kk=0; kk<NUM_SEGS; ++kk) {
                if(mask[kk]) {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj)
                        proposed2.second[jj] = iter->first.first[jj];
                } else {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj)
                        proposed2.second[jj] = iter->first.second[jj];
                }
            }

            char reassort_detail2 = 0;
            for(uint ll=0; ll < mask.size(); ++ll) {
                if(mask[ll]) {
                    reassort_detail2 |= (0x1 << ll);
                }
            }
            assert(reassort_detail2 != 0);

            // add first proposed vertex
            if(graph.vtoid.find(proposed1) != graph.vtoid.end() && proposed1 != iter->first) {
                uint vid1 = graph.vtoid.at(iter->first);
                uint vid2 = graph.vtoid.at(proposed1);
                Level2_graph::vv_to_ids::const_iterator edge_iter;
                if((edge_iter = graph.vvtoids.find(pair<uint,uint>(vid1, vid2))) == graph.vvtoids.end()) {
                    graph.add_edge(vid1, vid2, Level2_edge_data(1, 1, pair<char,char>(0, reassort_detail1)));
                    ++num_reassort_edges;
                } else {
                    Level2_graph::edge_type const &edge_data =  graph.idtoe.at(edge_iter->second.front());
                    if(edge_data.type == 0) {
                        if(static_cast<double>(edge_data.weight) * graph.graph_data.edge_costs.at(0) > graph.graph_data.edge_costs.at(1)) {
                            // if mutation is more expensive
                            // replace with reassortment
                            graph.remove_edge(edge_iter->second.front());
                            graph.add_edge(vid1, vid2, Level2_edge_data(1, 1, pair<char,char>(0, reassort_detail1)));
                            ++num_reassort_edges;
                        }
                    } else {
                        assert(edge_data.type == 1);
                    }
                }
            }

            // add second proposed vertex
            if(graph.vtoid.find(proposed2) != graph.vtoid.end() && proposed2 != iter->first) {
                uint vid1 = graph.vtoid.at(iter->first);
                uint vid2 = graph.vtoid.at(proposed2);
                Level2_graph::vv_to_ids::const_iterator edge_iter;
                if((edge_iter = graph.vvtoids.find(pair<uint,uint>(vid1, vid2))) == graph.vvtoids.end()) {
                    graph.add_edge(vid1, vid2, Level2_edge_data(1, 1, pair<char,char>(1, reassort_detail2)));
                    ++num_reassort_edges;
                } else {
                    Level2_graph::edge_type const &edge_data =  graph.idtoe.at(edge_iter->second.front());
                    if(edge_data.type == 0) {
                        if(static_cast<double>(edge_data.weight) * graph.graph_data.edge_costs.at(0) > graph.graph_data.edge_costs.at(1)) {
                            // if mutation is more expensive
                            graph.remove_edge(edge_iter->second.front());
                            graph.add_edge(vid1, vid2, Level2_edge_data(1, 1, pair<char,char>(1, reassort_detail2)));
                            ++num_reassort_edges;
                        }
                    } else {
                        assert(edge_data.type == 1);
                    }
                }
            }
        } // end while
    }
    cout << "Number of reassortment edges: " << num_reassort_edges << endl;
    assert(graph.check_multiedge() == 0);

    return 0;
}

//int Buneman::get_subset(bvector& mask) {
//    uint ii;
//    for(ii=0; ii<=mask.size(); ++ii) {
//        if(mask[ii]) {
//            mask[ii] = false;
//        } else {
//            mask[ii] = true;
//            break;
//        }
//    }
//    if(ii == mask.size()) {
//        return 0; 
//    }
//    return 1;
//}

template<typename VertexT, typename EdgeT, typename GraphDataT> 
int Buneman::dfs_search(Graph<VertexT, EdgeT, GraphDataT> const& graph, VertexT const& vert, boost::unordered_set<VertexT> & v_set, int& num_vert_reached) {
    if(v_set.find(vert) == v_set.end()) {
        v_set.insert(vert);
        ++num_vert_reached;
        typename Graph<VertexT, EdgeT, GraphDataT>::edge_list const &outedges = graph.outedges(vert);
        //cout << "num outgoing edges: " << outedges.size() << endl;
        for(typename Graph<VertexT, EdgeT, GraphDataT>::edge_list::const_iterator eiter = outedges.begin(); eiter != outedges.end(); ++eiter) {
            dfs_search(graph, graph.idtov.at(eiter->first), v_set, num_vert_reached); 
        }
    }
    return 0;
}

template<class T>
int Buneman::create_reassort_ILP(Graph<T>& graph) {
    // create MIP instance for reassortment ILP 
    string lp_name = "Level2_ILP";
    cout << lp_name << endl;
    glp_prob *lp = glp_create_prob();
    glp_set_prob_name(lp, lp_name.c_str());

    //////////////////////////
    // add variables of lp
    //////////////////////////

    // add s variables for every edge
    boost::unordered_map<boost::tuple<uint,uint,char>, int> map_edge_svar;
    boost::unordered_map<int, boost::tuple<int,int>  > map_edge_type;
    //vector<boost::unordered_map<pair<uint,uint>, int> > map_edge_fvar;
    //map_edge_fvar.resize(main_graph.terminal_vertices.size());

    for(uint v_id=0; v_id < graph.adj.size(); ++v_id) {
        for(edge_list::const_iterator eiter = graph.adj[v_id].first.begin(); eiter != graph.adj[v_id].first.end(); ++eiter) {
            int s_id = glp_add_cols(lp, 1);
            glp_set_col_kind(lp, s_id, GLP_BV); // binary
            //glp_set_col_bnds(lp, s_id, GLP_DB, 0.0, 1.0); // debug bounds
            map_edge_svar[boost::tuple<uint,uint,char>(v_id,eiter->get<0>(),eiter->get<1>())] = s_id;
            //map_edge_type[pair<uint,uint>(v_id,eiter->first)] = eiter->second;

            double cost;
            // set objective function
            if(eiter->get<1>() == 2 || eiter->get<1>() == 3) {
                // reassortment cost
                cost = graph.edge_costs[eiter->get<1>()];
                glp_set_obj_coef(lp, s_id, cost);
            } else {
                // mutation cost
                cost = (double)eiter->get<3>() * graph.edge_costs[eiter->get<1>()];
                glp_set_obj_coef(lp, s_id, cost);
            }
            map_edge_type[s_id] = boost::tuple<int,int>(eiter->get<2>(), eiter->get<3>());
        }
    }

    // arbitrarily choose the first terminal vertex to be the root
    assert(graph.vtoid.size() > 1);
    uint root_id = graph.vtoid[graph.terminal_vertices[0]];
    cout << "root id: " << root_id << endl;

    // For flow variables for terminal t
    // Don't include root in this
    for(uint term_vert=1; term_vert < graph.terminal_vertices.size(); ++term_vert) {
        cout << "processing terminal vertex: " << graph.terminal_vertices[term_vert] << ", " << graph.vtoid[graph.terminal_vertices[term_vert]] << endl;
        boost::unordered_map<boost::tuple<uint, uint, char>, int> map_edge_flow;
        // vert_id is the vertex id of the terminal vertex
        uint vert_id = graph.vtoid[graph.terminal_vertices[term_vert]];

        // create flow variables for term_vert
        for(uint v_id=0; v_id < graph.adj.size(); ++v_id) {
            // iterate over all edges in graph
            for(edge_list::const_iterator eiter = graph.adj[v_id].first.begin(); eiter != graph.adj[v_id].first.end(); ++eiter) {
                assert(map_edge_svar.find(boost::tuple<uint,uint,char>(v_id, eiter->get<0>(), eiter->get<1>())) != map_edge_svar.end());
                int s_id = map_edge_svar[boost::tuple<uint,uint,char>(v_id, eiter->get<0>(), eiter->get<1>())];
                int flow_id = glp_add_cols(lp, 1);
                //map_edge_fvar[term_vert][pair<uint,uint>(v_id, eiter->first)] = flow_id;
                map_edge_flow[boost::tuple<uint, uint, char>(v_id, eiter->get<0>(), eiter->get<1>())] = flow_id;
                // add column constraint: 0 <= f_uv
                // only one of the below is necessary, either binary or non-negative
                // make flow binary
                glp_set_col_kind(lp, flow_id, GLP_BV);
                //glp_set_col_bnds(lp, flow_id, GLP_LO, 0.0, 0.0);

                // add row constraint: f_uv <= s_uv
                //                     0 <= s_uv - f_uv
                int row_id = glp_add_rows(lp, 1);
                glp_set_row_bnds(lp, row_id, GLP_LO, 0.0, 0.0);
                int ind[3] = {0, s_id, flow_id};
                double val[3] = {0.0, 1.0, -1.0};
                glp_set_mat_row(lp, row_id, 2, ind, val);
            }
        } // done adding flow variables and s constraints for flow t

        // add flow conservation constraints for flow (t, extra_flow)
        for(uint v_id=0; v_id < graph.adj.size(); ++v_id) {
            // don't need flow conservation for root
            if(v_id == root_id) continue;
            // don't need flow conservation for terminal t
            if(v_id == vert_id) continue;

            // for v_id, create flow constraints for vertex
            boost::unordered_set<int> debug_set;
            int conserve_id = glp_add_rows(lp, 1);
            glp_set_row_bnds(lp, conserve_id, GLP_FX, 0.0, 0.0);

            // number of constraints is number of outgoing edges
            //    + incoming edges
            int len = graph.adj[v_id].first.size() + graph.adj[v_id].second.size();
            int ind[len+1];
            double val[len+1];

            // iterate over outgoing edges
            int index_counter = 1;
            for(edge_list::const_iterator eiter = graph.adj[v_id].first.begin(); eiter != graph.adj[v_id].first.end(); ++eiter) {
                uint fromedge = v_id;
                uint toedge = eiter->get<0>();
                // outgoing edge
                assert(map_edge_flow.find(boost::tuple<uint,uint,char>(fromedge,toedge, eiter->get<1>())) != map_edge_flow.end());
                ind[index_counter] = map_edge_flow[boost::tuple<uint,uint,int>(fromedge, toedge, eiter->get<1>())];
                debug_set.insert(ind[index_counter]);
                val[index_counter] = 1.0;
                ++index_counter;
            }
            for(edge_list::const_iterator eiter = graph.adj[v_id].second.begin(); eiter != graph.adj[v_id].second.end(); ++eiter) {
                uint fromedge = eiter->get<0>(); 
                uint toedge = v_id; 
                // incoming edge
                assert(map_edge_flow.find(boost::tuple<uint,uint,char>(fromedge,toedge, eiter->get<1>())) != map_edge_flow.end());
                ind[index_counter] = map_edge_flow[boost::tuple<uint,uint,char>(fromedge, toedge, eiter->get<1>())];
                if(debug_set.find(ind[index_counter]) != debug_set.end()) {
                    cout << "duplicate index" << endl;
                    cout << "edge: [" << graph.idtov[fromedge] << ", " << graph.idtov[toedge] << "]" << endl;
                    cout << "    : [" << fromedge << ", " << toedge << "]" << endl;
                    assert(false);
                }
                debug_set.insert(ind[index_counter]);
                val[index_counter] = -1.0;
                ++index_counter;
            }
            assert(index_counter == len+1);

            glp_set_mat_row(lp, conserve_id, len, ind, val);

        } // done adding flow conservation constraints for flow for term t

        //cout << "terminal constraints" << endl;
        // add constraint to prevent flow from leaving terminal t
        int term_id = glp_add_rows(lp, 1);
        glp_set_row_bnds(lp, term_id, GLP_FX, 0.0, 0.0);

        int len = graph.adj[vert_id].first.size();

        int ind[len+1]; 
        double val[len+1];

        int index_counter = 1;
        // go over all outgoing edges from terminal t
        for(edge_list::const_iterator eiter = graph.adj[vert_id].first.begin(); eiter != graph.adj[vert_id].first.end(); ++eiter) {
            uint fromedge = vert_id;
            uint toedge = eiter->get<0>();
            // outgoing edge
            assert(map_edge_flow.find(boost::tuple<uint,uint,char>(fromedge,toedge, eiter->get<1>())) != map_edge_flow.end());
            ind[index_counter] = map_edge_flow[boost::tuple<uint,uint,char>(fromedge, toedge, eiter->get<1>())];
            val[index_counter] = 1.0;
            ++index_counter;
        }
        assert(index_counter == len+1);

        glp_set_mat_row(lp, term_id, len, ind, val);
        // done adding term outflow constraints for flow t 

        // add constraint so that terminal receives 1 unit of flow 
        int term_in_id = glp_add_rows(lp, 1);
        glp_set_row_bnds(lp, term_in_id, GLP_FX, 1.0, 1.0);

        int len_in = graph.adj[vert_id].second.size();

        int ind_in[len_in+1];
        double val_in[len_in+1];

        index_counter = 1;
        // go over all incoming edges into terminal t 
        for(edge_list::const_iterator eiter = graph.adj[vert_id].second.begin(); eiter != graph.adj[vert_id].second.end(); ++eiter) {
            uint fromedge = eiter->get<0>();
            uint toedge = vert_id; 
            // incoming edge
            assert(map_edge_flow.find(boost::tuple<uint,uint,char>(fromedge,toedge, eiter->get<1>())) != map_edge_flow.end());
            ind_in[index_counter] = map_edge_flow[boost::tuple<uint,uint,char>(fromedge, toedge, eiter->get<1>())];
            val_in[index_counter] = 1.0;
            ++index_counter;
        }
        assert(index_counter == len_in+1);

        glp_set_mat_row(lp, term_in_id, len_in, ind_in, val_in);
        // done adding terminal inflow constraints

    } // done processing the flow for terminal term_vert

    // set parameters to solve solve mip
    glp_iocp parm; glp_init_iocp(&parm); // init parameters struct 
    parm.presolve = GLP_ON; // turn on presolver
    parm.msg_lev = GLP_MSG_ALL; // output verbosity

    assert(glp_write_lp(lp, NULL, "reassort.lp")==0);

    // write out mapping between svariable ids and corresponding edges
    string map_file_name = "svarmapping.map";
    cout << "Writing map data to `" << map_file_name << "\'..." << endl;
    ofstream fh(map_file_name.c_str());
    int lines_written = 0;
    // write out root
    fh << graph.idtov[root_id] << '\n';
    ++lines_written;

    assert(graph.terminal_vertices.size() > 1);
    // write out terminals
    for(uint ii=1; ii < graph.terminal_vertices.size() - 1; ++ii) {
        fh << graph.terminal_vertices[ii] << '@'; 
    }
    fh << graph.terminal_vertices[graph.terminal_vertices.size()-1];
    fh << '\n';
    ++lines_written;

    // write out segment positions
    for(uint ii=0; ii<NUM_SEGS-1; ++ii) {
        fh << segment_start[ii] << ' ' << segment_end[ii] << '@';
    }
    fh << segment_start[NUM_SEGS-1] << ' ' << segment_end[NUM_SEGS-1] << '\n';

    for(boost::unordered_map<boost::tuple<uint,uint,char>,int>::const_iterator siter=map_edge_svar.begin(); siter!=map_edge_svar.end(); ++siter) {
        //if(map_edge_type[siter->first] == -2 || map_edge_type[siter->first] >= 0 ) {
        fh << siter->second << '@' << graph.idtov[siter->first.get<0>()] << '@' << graph.idtov[siter->first.get<1>()] << '@' << (int)siter->first.get<2>() << '@' << map_edge_type[siter->second].get<0>() << '@' << map_edge_type[siter->second].get<1>() << '\n';
        ++lines_written;
        //}
    }
    cout << lines_written << " lines were written" << endl;
    fh.close();

    //cout << "cur_graph.num_vertices(): " << cur_graph.num_vertices() << endl;
    //cout << "cur_graph.terminal_vertices.size(): " << cur_graph.terminal_vertices.size() << endl;

    // Use glpk to solve the ILP
    //cout << "Solving ILP" << endl;
    //int ret = glp_intopt(lp, &parm);
    //assert(ret == 0);
    //cout << "ret: " << ret << endl;
    //cout << "obj value: " << glp_mip_obj_val(lp) << endl;

    // display solution, s vars
    //cout << "displaying solution:" << endl;
    //cout << "svars" << endl;
    //for(boost::unordered_map<pair<uint,uint>, int>::const_iterator siter = map_edge_svar.begin(); siter != map_edge_svar.end(); ++siter) {
    //    double sol_value = glp_mip_col_val(lp, siter->second);
    //    if(sol_value != 0) {
    //        cout << "edge: ";
    //        cout << main_graph.idtov[siter->first.first] << " " << main_graph.idtov[siter->first.second] << endl;
    //        cout << "val:  " << sol_value << endl;
    //    }
    //}

    // display solution, f vars
    //cout << "flow vars" << endl;
    //cout << "root: " << cur_graph.terminal_vertices[0] << endl;
    //for(uint tt=1; tt<cur_graph.terminal_vertices.size(); ++tt) {
    //    cout << "flow for terminal " << tt << endl;
    //    cout << "terminal: " << cur_graph.terminal_vertices[tt] << endl;
    //    for(boost::unordered_map<pair<uint,uint>, int>::const_iterator fiter = map_edge_fvar[tt].begin(); fiter != map_edge_fvar[tt].end(); ++fiter) {
    //        double sol_value = glp_mip_col_val(lp, fiter->second);
    //        if(sol_value != 0) {
    //            cout << "edge: ";
    //            cout << cur_graph.idtov[fiter->first.first] << " " << cur_graph.idtov[fiter->first.second] << endl;
    //            cout << "val:  " << sol_value << endl;
    //        }
    //    }
    //}

    //int obj_value = glp_mip_obj_val(lp);
    //cout << "obj_value: " << obj_value << endl;
    // this is a bad assertion since steiner tree might use
    // intermediate vertices
    //assert(obj_value == static_cast<int>(cur_graph.terminal_vertices.size()) - 1);


    // delete lp
    glp_delete_prob(lp);
    // done with ILP with reassortment vertices
    return 0;
}

int Buneman::create_big_graph() {
    cout << "adding vertices" << endl;
    // add vertices to big graph
    for(Graph<cvector>::v_to_id::const_iterator iter1 = main_graph.vtoid.begin(); iter1 != main_graph.vtoid.end(); ++iter1) {
        for(Graph<cvector>::v_to_id::const_iterator iter2 = main_graph.vtoid.begin(); iter2 != main_graph.vtoid.end(); ++iter2) {
            pair<uint, bool> add_ret = big_graph.add_unique_vertex(pair<cvector,cvector>(iter1->first, iter2->first));
            assert(add_ret.second); // assert that a unique vertex was added
        }
    }
    assert(big_graph.num_vertices() == main_graph.num_vertices() * main_graph.num_vertices());
    cout << "adding edges" << endl;

    // add edges to big graph
    // mutate first hap 
    uint num_sites = big_graph.vtoid.begin()->first.first.size();
    for(Graph<pair<cvector,cvector> >::v_to_id::const_iterator iter1 = big_graph.vtoid.begin(); iter1 != big_graph.vtoid.end(); ++iter1) {
        for(uint mm=0; mm < num_sites; ++mm) {
            pair<cvector, cvector> proposed = iter1->first;
            proposed.first[mm] = 1 - proposed.first[mm]; 

            if(big_graph.vtoid.find(proposed) != big_graph.vtoid.end()) {
                // edge type should be different in final version
                bool add_ret = big_graph.add_edge(iter1->first, proposed, 0, mm, 1);
                assert(add_ret);
            }
        }
    }

    // mutate second hap 
    assert(big_graph.vtoid.begin()->first.first.size() == big_graph.vtoid.begin()->first.second.size());
    for(Graph<pair<cvector,cvector> >::v_to_id::const_iterator iter2 = big_graph.vtoid.begin(); iter2 != big_graph.vtoid.end(); ++iter2) {
        for(uint mm=0; mm < num_sites; ++mm) {
            pair<cvector, cvector> proposed = iter2->first;
            proposed.second[mm] = 1 - proposed.second[mm]; 

            if(big_graph.vtoid.find(proposed) != big_graph.vtoid.end()) {
                // edge type should be different in final version
                bool add_ret = big_graph.add_edge(iter2->first, proposed, 1, mm, 1);
                assert(add_ret);
            }
        }
    }

    /////////////////////////
    // reassortment edges
    for(Graph<pair<cvector,cvector> >::v_to_id::const_iterator iter = big_graph.vtoid.begin(); iter != big_graph.vtoid.end(); ++iter) {
        vector<bool> mask(NUM_SEGS, false);
        while(get_subset(mask)) {
            bool allones = true;
            for(uint ii=0; ii < mask.size(); ++ii) {
                if(mask[ii] == false) {
                    allones = false;
                    break;
                }
            }
            if(allones) continue;
            // reassort first hap
            pair<cvector, cvector> proposed = iter->first;
            for(uint kk=0; kk<NUM_SEGS; ++kk) {
                if(mask[kk]) {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj) {
                        proposed.first[jj] = iter->first.first[jj];
                    }
                } else {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj) {
                        proposed.first[jj] = iter->first.second[jj];
                    }
                }
            }

            int reassort_detail = 0;
            for(uint ll=0; ll < mask.size(); ++ll) {
                if(mask[ll]) {
                    reassort_detail |= (0x1 << ll);
                }
            }
            assert(reassort_detail != 0);

            if(big_graph.vtoid.find(proposed) != big_graph.vtoid.end() && proposed != iter->first) {
                big_graph.add_edge(iter->first, proposed, 2, reassort_detail, 1);
            }

            // reassort second hap
            pair<cvector, cvector> proposed2 = iter->first;
            for(uint kk=0; kk<NUM_SEGS; ++kk) {
                if(mask[kk]) {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj) {
                        proposed2.second[jj] = iter->first.first[jj];
                    }
                } else {
                    for(int jj=segment_start[kk]; jj<segment_end[kk]; ++jj) {
                        proposed2.second[jj] = iter->first.second[jj];
                    }
                }
            }

            reassort_detail = 0;
            for(uint ll=0; ll < mask.size(); ++ll) {
                if(mask[ll]) {
                    reassort_detail |= (0x1 << ll);
                }
            }
            assert(reassort_detail != 0);

            if(big_graph.vtoid.find(proposed2) != big_graph.vtoid.end() && proposed2 != iter->first) {
                big_graph.add_edge(iter->first, proposed2, 3, reassort_detail, 1);
            }
        }

        // add in special reassortment where the child agrees
        //    completely with one of the parents
        // the edge cost is 0 in this case
        // only do this if the parents are different from each other
        if(iter->first.first != iter->first.second) {
            int reassort_detail = 0;
            pair<cvector, cvector> childleft = pair<cvector, cvector>(iter->first.first, iter->first.first);
            big_graph.add_edge(iter->first, childleft, 4, reassort_detail, 1);

            pair<cvector, cvector> childright = pair<cvector, cvector>(iter->first.second, iter->first.second);
            big_graph.add_edge(iter->first, childright, 5, reassort_detail, 1);
        }
    }

    // create list of terminal vertices
    big_graph.terminal_vertices.clear();
    for(vector<cvector>::const_iterator viter=main_graph.terminal_vertices.begin(); viter != main_graph.terminal_vertices.end(); ++viter) {
        big_graph.terminal_vertices.push_back(pair<cvector,cvector>(*viter, *viter));

        assert(big_graph.terminal_set.find(pair<cvector,cvector>(*viter,*viter)) == big_graph.terminal_set.end());
        big_graph.terminal_set.insert(pair<cvector,cvector>(*viter, *viter));
    }

    return 0;
}

namespace boost {
    std::size_t hash_value(boost::tuple<uint,uint,int> const& t) {
        std::size_t seed = 0;
        boost::hash_combine(seed, boost::get<0>(t));
        boost::hash_combine(seed, boost::get<1>(t));
        boost::hash_combine(seed, boost::get<2>(t));
        return seed;
    }

    std::size_t hash_value(boost::tuple<uint,uint,char> const& t) {
        std::size_t seed = 0;
        boost::hash_combine(seed, boost::get<0>(t));
        boost::hash_combine(seed, boost::get<1>(t));
        boost::hash_combine(seed, boost::get<2>(t));
        return seed;
    }

    std::size_t hash_value(boost::tuple<pair<cvector,cvector>, pair<cvector,cvector>, char> const& t) {
        std::size_t seed = 0;
        boost::hash_combine(seed, boost::get<0>(t));
        boost::hash_combine(seed, boost::get<1>(t));
        boost::hash_combine(seed, boost::get<2>(t));
        return seed;
    }

    std::size_t hash_value(boost::tuple<cvector, cvector, char> const& t) {
        std::size_t seed = 0;
        boost::hash_combine(seed, boost::get<0>(t));
        boost::hash_combine(seed, boost::get<1>(t));
        boost::hash_combine(seed, boost::get<2>(t));
        return seed;
    }
}
